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Photosystem II (PSII) and its associated light-harvesting complex II (LHCII) are highly concen- 
trated in the stacked grana regions of photosynthetic thylakoid membranes. Within the membrane, 
PSII-LHCII supercomplexes can be arranged in disordered packings, ordered arrays, or mixtures 
thereof. The physical driving forces underlying array formation are unknown, complicating attempts 
to determine a possible functional role for arrays in regulating light harvesting or energy conversion 
efficiency. Here we introduce a coarse-grained model of protein interactions in coupled photosyn- 
thetic membranes, focusing on just two particle types that feature simple shapes and potential 
energies motivated by structural studies. Reporting on computer simulations of the model's equi- 
librium fluctuations, we demonstrate its success in reproducing diverse structural features observed 
in experiments, including extended PSII-LHCII arrays. Free energy calculations reveal that the 
appearance of arrays marks a phase transition from the disordered fluid state to a system-spanning 
crystal, which can easily be arrested by thermodynamic constraints or slow dynamics. The region 
of fluid-crystal coexistence is broad, encompassing much of the physiologically relevant parameter 
regime. Our results suggest that grana membranes lie at or near phase coexistence, conferring 
significant structural and functional flexibility to this densely packed membrane protein system. 



I. INTRODUCTION 

Photosynthetic efficiency relies on precise spatial or- 
ganization of sites for light harvesting, exciton trans- 
port, and charge separation 1 . In higher plants and 
green algae, these functions are carried out by pigment- 
proteins in the thylakoid membrane, with PSII host- 
ing charge separation and oxygen evolution functionality 
and LHCII acting as its associated light-harvesting an- 
tenna. PSII and LHCII exhibit self- assembled organiza- 
tion on a range of length scales [21 E]. At the macro- 
molecular scale, dimers of PSII core complexes asso- 
ciate specifically with 1-6 trimers of LHCII and up to 
2 monomers of each of the minor light harvesting com- 
plexes (CP24, CP26, and CP29) to form a family of PSII 
supercomplexes [1HS]. On larger scales, the thylakoid 
membrane is differentiated into stacked discs of tightly 
appressed membrane called grana, tubes or sheets of un- 
appressed membrane called stroma lamellae, and contro- 
versial junctional regions called margins [5J[7HS]- PSII 
and LHCII are among the many proteins that display 
dramatic lateral heterogeneity within this complex mem- 
brane architecture — PSII and LHCII typically localize 
to appressed grana membranes, photosystem I (with its 
light-harvesting complex I) and ATP synthase are pre- 
dominantly found in unappressed membranes, and other 
proteins such as cytochrome b e f are present in both re- 
gions JTUl HI] ■ The degrees of both supercomplex forma- 
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tion and lateral heterogeneity are dynamically regulated 
in reponse to environmental factors by the processes of 
state transitions and PSII repair, in which phosphoryla- 
tion of LHCII or PSII (respectively) is correlated with 
partial or complete dissolution of the supercomplex and 
migration of the phosphorylated protein species toward 
stroma lamellae or margin regions; many mechanistic de- 
tails remain under debate [L2TTT5] . 

Within grana stacks, the most striking features of pro- 
tein organization are regular 2D patterns termed arrays. 
Since their early description in the 1960s |16j . these 
motifs have frequently beeen observed in electron mi- 
croscopy and atomic force microscopy studies of stacked 
thylakoid membranes, and are composed of tens or hun- 
dreds of unit cells containing one PSII core dimer and 
variable quantities of LHCII (reviewed in [21 [TTl [18]). 
Array extent and unit cell have been correlated with a 
number of experimental factors, including mutations of 
light-harvesting complex [T9H23] , PsbS [24] , or photosys- 
tem I [25] proteins; cold storage 26, 27]; and acclimation 
to low light [3SJ Yet, in part because observed ar- 
rays can vary widely between similar samples [501 15T] . 
quantitatively predictive statements about the structure 
and function of PSII arrays have not emerged. Hypothe- 
sized functional rationales for PSII arrays primarily focus 
on predicted exciton transport effects, such as optimiz- 
ing antenna size, excitation quenching site availability, 
or reaction center connectivity to match light conditions 
[521 [55] , and on protein and electron carrier mobility ar- 
guments [5H |5S] , but no functional role has been proven. 

Nanoscale maps of the relative arrangement of LHCII 
and PSII in native-like grana membrane environments 
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will be necessary for building a complete picture of exci- 
ton flow during photosynthetic light harvesting and en- 
ergy conversion. The oxygen-evolving complex of PSII 
protrudes above the plane of the membrane on the lu- 
menal side [5t)H10] . allowing it to be easily identified in 
electron microscopy and atomic force microscopy images 
of grana membranes. LHCII lacks a lumenal protrusion 
[32 [32] ; this has prevented microscopists from simulta- 
neously mapping the positions of LHCII and PSII in dis- 
ordered membrane environments (i.e., outside of arrays). 
Super-resolution fluorescence imaging techniques are a 
promising alternative, but have not yet been successfully 
applied to grana membranes, which are already highly 
fluorescent. Computer simulation has the potential to 
complement experiment and address this important gap 
in structural and mechanistic understanding. 

The first computational study of thylakoid protein or- 
ganization focused on a single membrane layer, modeled 
as a flat 2D sheet, and treated grana and stroma lamellae 
membrane proteins as hard disc-shaped particles moving 
within that sheet [33] . Tremmel and coworkers developed 
a distinct model restricted to grana membrane protein 
complexes (LHCII trimers, PSII supercomplexes, and 
cytochrome b 6 f) that used single-particle Monte Carlo 
moves to sample configurations of hard particles with 
highly detailed shapes on a single 2D lattice [331 US]; 
they later extended the model by introducing short-range 
patchy interactions between protein particles |46j , a stan- 
dard means of modeling anisotropic protein-protein inter- 
actions [7fH5Tj . None of these publications reported the 
emergence of PSII arrays in computer simulation. 

In this study, we use a simple nanoscale computational 
model of LHCII and PSII in stacked membranes to an- 
alyze grana protein organization. In particular, we ex- 
amined PSII array formation in system sizes comparable 
to full grana discs (hundreds of nanometers in diameter) . 
We extend previous computational approaches in two key 
ways: we simulate multiple, coupled membrane layers in 
a grana stack; and we employ computational methods 
capable of exploring highly cooperative transitions. By 
thoroughly examining thermal fluctuations in our model 
of a PSII-LHCII binary mixture, we reveal a phase tran- 
sition between a relatively dilute, disordered PSII-LHCII 
fluid and a dense, ordered PSII-LHCII crystal. Physio- 
logical protein concentrations are at fluid-crystal phase 
coexistence or near the boundary between the fluid and 
coexistence regions, where small changes in protein den- 
sity or interactions can lead to dramatic shifts in the 
observed degree of array formation and in any array- 
dependent functionality. 



II. RESULTS 

A. Model is founded on in vivo phenomenology 

Model details are illustrated in Fig. [T] and Fig. SI, 
and further described in the Methods and Text SI. 




LHCII-LHCII interlayer (stacking) attraction 




FIG. 1: Graphical summary of the model, a: Side 
view cartoon of protein geometry and interactions. PSII 
(blue) and LHCII (green) particles reside in two coupled lay- 
ers, each of which represents a lipid bilayer grana membrane 
(tan) surrounded by stromal and lumenal aqueous regions 
(pale blue), b and c: Top view from the lumenal side of 
the coarse-graining procedure. Solved protein structures of 
LHCII trimers are roughly circular; model LHCII particles 
are hard discs with diameter 6.5 nm (dark green: S-LHCII; 
medium green: M-LHCII). PSII C2S2 supercomplexes are 
roughly rod-shaped; model PSII particles are hard discorect- 
angles with diameter 12 nm and tip-to-tip length 26.5 nm 
(blue). The orange dots show LHCII binding sites that allow 
formation of supramolecular complexes, including (0282)2 (b, 
governed by esl-p) and C2S2M2 (c, governed by eml-p). 
Protein structures adapted from [2] [5]. d: Potential energy of 
the LHCII interlayer stacking interaction. Insets: side view 
cartoons of unbound LHCII (light green) in stacked (left) and 
unstacked (right) configurations. 

Our model consists of two coupled 2D layers, represent- 
ing stroma-paired grana membranes, which comprise the 
minimal system for PSII array formation [52] . These two 
layers can be thought of as an isolated membrane pair, or 
as two stroma-paired membranes in a larger grana stack 
where lumen-side interactions are assumed to be negligi- 
ble [2J. Two types of hard particles reside in these lay- 
ers: disc-shaped species that represent LHCII trimers, 
and rod-shaped species that represent PSII supercom- 
plexes (specifically the so-called C2S2 supercomplex in 
Chlamydomonas [5], spinach |53| . and Arabidopsis [S]). 
PSII particles, like PSII protein supercomplexes, have 
two constituent embedded LHCIIs that move with the 
PSII particle as a rigid body (Fig. [T|d and c, Methods). 
The embedded LHCIIs are constrained to their physiolog- 
ical locations on either side of the PSII axis [S3], making 
the PSII particles chiral (Fig. SI). 

In addition to excluding area within the layer that they 
occupy, these particles attract one another in three dis- 
tinct ways. One mode of attraction acts between LHCII 
particles in opposing layers. Electrostatic interactions 
between LHCII stromal faces contribute to grana stack 
integrity |55j and are implicated in PSII array formation 
in experiments [301 152). Standfuss and coworkers have 
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proposed a "velcro-like" qualitative model for the LHCII 
stacking interaction [42] , but a quantitative understand- 
ing is lacking. We have chosen a simple phcnomcnological 
form for the LHCII-LHCII interlayer attraction that fa- 
vors LHCII face-to-face contact (Fig. [lp,d). We set the 
strength of the attraction to £l-l = 4 feaT, which cre- 
ates strong but reversible binding and is consistent with 
electrostatic measurements (see Text SI). All LHCII par- 
ticles can participate in this interaction. 



The other attractions associate PSII and LHCII parti- 
cles residing in the same layer. When protein complexes 
are isolated from thylakoid membranes of Arabidopsis, 
spinach, or Chlamydomonas, LHCII complexes are found 
to be distributed among at least three different local en- 
vironments: they can be strongly bound within a su- 
percomplex (embedded or S-LHCII), moderately bound 
to the edge of a supercomplex (bound or M-LHCII), or 
structurally detached from PSII (free LHCII) [1 G3 H EJ • 
Motivated by the apparent equilibrium between moder- 
ately bound and free LHCIIs, we introduce two interac- 
tion sites on the periphery of each PSII particle, as illus- 
trated by the orange dots in Fig. [T]d and c and Fig. SI. 
These interaction sites attract LHCIIs in the same layer 
over a short range (~1 nm), such that modeled M-LHCIIs 
in the "bound" state are constrained to locations relative 
to PSII that are consistent with published C2S2M.,; com- 
plexes [HE]. 

Experiments suggest that the intralayer interactions 
between PSII supercomplexes and free LHCIIs are simi- 
lar to yet distinct from those between a PSII supercom- 
plex and the S-LHCII of another supercomplex; for in- 
stance, binding of M-LHCII (Fig. [lj;) appears to require 
PSII subunit CP24 in Arabidopsis, while association of 
S-LHCII with the same site (Fig. ^p) may not [21]. We 
therefore define separate energy scales for these interac- 
tions (cml-p and £sl-p, respectively). Based on mea- 
surements of intramembrane associations between other 
protein complexes in photosynthetic membranes |56j . 
these energies are expected to be modest, on the order of 
one or a few k^T; we focus on £sl-p = £ml-p = Zk^,T. 

We emphasize that this model includes only the short- 
ranged intramembrane protein-protein attractions that 
are most conclusively established by single-particle struc- 
tural studies. Longer-range patterns in our simulations 
can only arise from emergent correlations involving many 
particles. We also emphasize that the values of the inter- 
action strengths that we have selected are intended to be 
representative of the physiology of a generic grana mem- 
brane, and not to represent any specific species or growth 
condition. Experimental measurements of the interaction 
strengths, which would be difficult but not inconceivable 
(e.g., along the lines of Ref. [56]). could support future 
efforts to tailor the model to precisely match specific ex- 
perimental conditions. 



B. Simulations capture experimentally-determined 
intramembrane protein organization 

We performed Monte Carlo simulations to sample equi- 
librium particle configurations at a wide range of particle 
concentrations within experimentally relevant parameter 
regimes. Two standard experimental measures of grana 
protein content are the LHCILPSII ratio (specifically, the 
mole ratio (j> of free LHCII trimers to PSII C2S2 super- 
complexes) and the number density of PSII (typically 
PSII per square micron of membrane). These metrics 
can be combined to estimate the total protein packing 
fraction p (i.e., the fraction of grana membrane area oc- 
cupied by proteins) . Published values of these quantities 
vary significantly based on plant growth conditions and 
membrane preparation protocols: 4> is typically in the 
range ~ 2-6 [251 [571 [5HJ , and p is typically in the range 
« 0.6-0.8 [291 SUES]. 
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FIG. 2: Intramembrane structure of a coupled pair of 
model grana membranes at 1750 PSII//im 2 . a: Snap- 
shot of the lumen-side-up membrane layer from a simulation 
at p — 0.70, (f> = 3.36 showing a representative disordered con- 
figuration. Example supramolecular complexes are circled, 
clockwise from top: C 2 S 2 , C 2 S 2 M, (C 2 S 2 ) 2 , C 2 S 2 M 2 . Color 
scheme as in Fig. [I] scale bar = 50 nm. b and c: Comparison 
between experimental and simulated data for the statistics 
of intramembrane PSII center-to-center separation distances, 
b shows the distribution of nearest-neighbor distances, and 
c shows the radial distribution function g(r). Freeze-fracture 
EM data on isolated spinach grana membranes [45] . with 1700 
PSII/pm 2 on average, are consistent with these simulations 
of the disordered state. In contrast, simulations at p — 0.75 
show PSII aggregation and ordering. 

Simulations of moderately dense systems (p = 0.70, 
4> = 3.36) showed good agreement with a variety of ex- 
perimental probes of intramembrane protein organiza- 
tion. Fig. [2] shows data from lumen-side-up layers of 
simulated membrane stacks, the perspective from which 
PSII positions are routinely detected experimentally [17] . 
LHCII and PSII particles were distributed uniformly and 
isotropically throughout the layer, with PSII orientations 
correlated over about 30 nm (2.5 PSII widths). Statis- 
tics of PSII-PSII separation distances matched well with 
experimental results, with computed and experimen- 
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tal nearest-neighbor PSII center-to-center distances of 
18.7±2.9 nm and 19.0±4.1 nm, respectively, and weakly 
structured radial distribution functions g(r) (Fig. and 
c, [IS]); small offsets in peak positions could be due to our 
simplified particle shapes. These correlations arise both 
from direct association between PSII pairs and from ef- 
fective forces mediated by LHCIIs and by particles in the 
opposing layer. Contributions from LHCII fluctuations 
include "depletion" attractions [5D], which tend to clus- 
ter PSIIs in order to maximize space available for the 
smaller LHCII species to explore. 

Though simple in form, our model of PSII-LHCII 
intralayer attraction was sufficient to stabilize sig- 
nificant populations of several known supramolecular 
complexes — PSII supercomplex dimers ((0282)2), PSII 
supercomplex with one additional LHCII (C2S2M), 
and PSII supercomplex with two additional LHCIIs 
(C2S2M2) — in equilibrium with a pool of lone PSII su- 
percomplexes (C2S2). The balance of this equilibrium 
can be tuned by changing the strength of the intralayer 
attraction; at e SL _ P = cml-p = 2 k B T, LHCII binding 
and unbinding is facile, and approximately 25% of LHCII 
particles are bound to PSII supercomplexes under these 
conditions. 



b = 3.36, p a 0.70 



4> = 4.22, p = 0.75 



LHCII stacking creates interlayer PSII 
correlations 



Our simple model of LHCII stacking interactions sim- 
ilarly proved sufficient to generate a range of struc- 
tural motifs that have been inferred from experiment. 
Specifically, we frequently observed (I) pairs of stacked, 
nearly parallel PSIIs, in which the LHCIIs chirally em- 
bedded in one PSII are both aligned with the embed- 
ded LHCIIs of the opposing PSII (of opposite chirality) 
(Fig. [3^); (2) nearly perpendicular pairs, in which all 
embedded LHCIIs stack over free LHCIIs in the other 
layer (Fig. |3Jd); and (3) rows of rotated PSIIs in one 
layer stacked atop a parallel but oppositely rotated row 
in the other layer, in which the two LHCIIs embedded in 
a given PSII are aligned with those of two different PSIIs 
in the opposing layer (Fig. |3j:) . In the absence of LHCII 
stacking interactions (i.e., when = 0), these PSII 

correlations disappeared (Fig. S2). 

All three of these motifs have been reported in sepa- 
rate experimental studies |30j [52] [53] ■ From those obser- 
vations, however, it was not clear whether different mo- 
tifs could be simultaneously abundant; nor could those 
authors establish LHCII stacking as a sufficient driving 
force for various interlayer correlations. Our results in- 
dicate that LHCII stacking can indeed drive all of these 
associations among PSIIs in opposing layers. 




FIG. 3: LHCII stacking effects on PSII organization. 

a-c: PSII interlayer motifs, as described in the text. Green 
and blue outlines: LHCII and PSII particles in the upper 
(lumen-side-up) layer; purple and red outlines: LHCII and 
PSII particles in the lower (stroma-side-up) layer. Black lines 
drawn parallel to the long axis of selected rods highlight orien- 
tational relationships between PSII particles in different lay- 
ers. Orientation correlations A9 also appear in Fig. S2c. Scale 
bars = 10 nm. d: Snapshot of the top layer from a simulation 
at <j> — 4.22, p = 0.75 showing a representative array. Color 
scheme as in Fig. T] except arrayed PSII are colored red (see 
Methods). Scale bar = 50 nm. e and f: Magnified views of 
the boxed region in panel d. Scale bars = 20 nm. The role 
of intralayer attractions is highlighted in e by showing only 
particles in the top layer and indicating the locations of PSII- 
LHCII interaction sites on each PSII. The stabilizing role of 
stacking is highlighted in f by showing particles from both 
layers in outline form (as in panels a-c). 



D. Simulated PSII arrays depend on packing 
fraction and attraction strength 



Above a packing fraction of p « 0.7, simulations ex- 
hibited sizable ordered arrays, featuring alternating rows 
of PSIIs and LHCIIs (Fig. [|l-f ) . Interestingly, PSIIs in 
these configurations do not engage in direct intralayer 
attractions. Instead, adhesion between PSII rows is pro- 
vided by intralayer attractions to the interspersed M- 
LHCIIs that bridge between PSII rows. Each row of 
PSIIs is stabilized by stacking interactions with a PSII 
row in the opposing layer, as in Fig. [3b (and by the less 
geometrically specific depletion attractions). This key 
role for stacking of embedded LHCIIs is highlighted by 
an absence of ordered arrays in simulations with £l-l = 
(Fig. S3). The experimental correlation between arrays 
and stacking [52] supports the realism of our coarse- 
grained model, as well as the conclusion that both modes 
of attraction in our model are essential for ordering. 
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E. A fluid-crystal phase transition is manifested in 
osmotic ensemble simulations 



Configurations in which tens or hundreds of PSIIs 
cluster tightly together, such as in Fig. [3] and in many 
EM images, reflect strong emergent forces of associa- 
tion, sufficient to offset the entropic cost of sequestering 
the constituent complexes from the surrounding disor- 
dered environment. Computer simulations allow us to 
address whether these large arrays further signify a more 
profound underlying phenomenon, namely, a true phase 
transition from a disordered "fluid" of relatively low PSII 
density to a system-spanning crystal of tightly packed 
PSIIs. The simulations we have described thus far are 
not suited to address this question, nor would be experi- 
ments probing isolated grana membranes with fixed pro- 
tein content. In these situations crystallization could not 
proceed to completion, simply because the system's net 
composition and total area are fixed at values inconsis- 
tent with the crystalline phase. Indeed, in simulations of 
closed systems like those shown in Fig. [3] array growth 
halted before all PSII had been incorporated, creating 
a dynamic equilibrium between arrayed and disordered 
PSII (Fig. S4). 

Phase transitions can be more readily identified by 
studying open systems such as the osmotic ensemble, in 
which area and composition can fluctuate subject to ex- 
ternal fields at fixed temperature 61 . First, we imposed 
a 2D "pressure" p that regulates changes in total area. 
More precisely, we fixed the osmotic pressure that par- 
ticles experience parallel to the plane of the membrane; 
simulated area fluctuations implicitly represent addition 
or subtraction of lipids from the membrane. Second, 
we allowed the population of one component (we chose 
LHCII) to fluctuate at fixed chemical potential p\,. Cor- 
responding number fluctuations in a real system involve 
exchanging material with a very large bath. In intact thy- 
lakoids, the stroma lamellae and other connected grana 
stacks could play the role of a bath. 

Our examination of detailed phase behavior focused 
on the type of array shown in Fig. [3]i by restricting 
the model interactions. In the simulations described be- 
low, intralayer attractions between PSIIs and embed- 
ded LHCIIs were omitted by setting £sl-p = 0, while 
attractions to free LHCIIs are retained by maintaining 
£ml-p = 2 fee? 1 ; thus, C2S2Ma; complexes were the only 
single-layer supramolecular structures directly stabilized 
by model energetics. Since the intralayer attractions 
we disabled do not directly contribute to the stability 
of these arrays, we expect the fully interacting model 
to exhibit qualitatively similar phase behavior. The 
prevalence of C2S2Ma; complexes in Arabidopsis [5] and 
Chlamydomonas [5] suggests that the restricted model 
may be particularly appropriate for these organisms. 

Varying pressure at fixed chemical potential /il pro- 
duced a sharp change in average packing fraction 
(Fig. |4^) , indicating a highly cooperative transition. The 
concomitantly sudden appearance of a system-spanning 



PSII array (Fig. [4;, Fig. S6) suggests that the degree of 
cooperativity wou d grow with system size, as in a first- 
order phase transition, and identifies the high-pressure 
phase as crystalline. For the finite, micron-scale sys- 
tem that we simulated, the jump in packing fraction p 
is necessarily rounded and could only become discontin- 
uous in the thermodynamic limit of an infinitely large 
system. Demonstrating true phase behavior would re- 
quire an analysis of scaling as this limit is approached. 
Given the limited spatial extent of natural thylakoids, 
we instead scrutinized remnant hallmarks of phase coex- 
istence in a finite system, specifically bistable free energy 
profiles and the presence of stable interfaces. 

We computed the free energy F(p) as a function of 
packing fraction for specific values of the thermodynamic 
parameters (p, fj,^, T) using umbrella sampling (Fig. S5). 
Free energy profiles at many other values of these param- 
eters were then calculated by thermodynamic reweight- 
ing [55]. These profiles exhibit two distinct basins, with 
minima at pf < 0.7 and p c > 0.8, over a range of 
pressures (Fig. [4Jd) . The low-packing fraction minimum 
corresponds to a disordered PSII-LHCII two-component 
fluid (bottom of Fig. |4p); configurations representative 
of the high-packing fraction minimum are nearly perfect 
PSII-LHCII co-crystals (top left of Fig.[4j Fig. S6). Con- 
straining the packing fraction to lie midway between pf 
and p c produces heterogeneous structures in which fluid 
and crystal coexist, separated by a system-spanning in- 
terface (Fig. [4J2) . Together with observations of hysteresis 
when pressure is cycled above and below the coexistence 
pressure p* (Fig. S7), these observations point strongly 
to a first-order phase transition. 



F. Phase coexistence region includes physiological 
conditions 

From free energy profiles like those of Fig. |4Jd, we con- 
structed a phase diagram in the plane of packing frac- 
tion p and composition For a given chemical potential 
the coexistence pressure p* can be uniquely identified 
as the pressure that maximizes the variance of packing 
fraction fluctuations; note the large standard deviation 
at p* in Fig. |4^. Values of p and <p f° r the two phases 
at the thermodynamic state (p^, p*) determine a pair of 
points at the boundaries of the crystal-fluid coexistence 
region. Repeating this procedure for different values of 
we traced out coexistence curves bounding regions 
where homogeneous fluid and crystal phases are thermo- 
dynamically stable (Fig. [5]). 

The dashed lines in Fig. [5] trace tie lines, i.e., lines 
along which the relative extent of the two pure phase 
regions varies while the packing fraction and composition 
of each phase remains constant (as determined by the 
cndpoints). These tie lines thus enable straightforward 
predictions for thermodynamic states in the coexistence 
region, requiring no characterization beyond the physical 
properties of the endpoints. 



G 




packing fraction p 



FIG. 4: Umbrella sampling simulations provide evidence for a phase transition, a: Applied pressure p versus 
packing fraction p along a line of constant chemical potential shows a sharp crossover at p* from a low-pressure, low-packing 
fraction regime to a high-pressure, high-packing fraction regime. Means (line and points) and root-mean-squared fluctuations 
(whiskers) of p at each pressure are calculated from probability distributions derived from free energy surfaces like those in 
panel b. Fluctuations, and therefore whiskers, are large in the vicinity of p*; for clarity we show only one whisker in this region, 
b: Free energy as a function of p for a system at relative chemical potential /2l =0.1 feT and three values of pressure near 
coexistence: within the fluid phase (blue, p — 0.06996 fceT/nm 2 ), within the crystalline phase (red, p = 0.07020 feT/nm 2 ), and 
at coexistence (green, p* = 0.07007 k^T/m-a?). Error bars estimated from the MBAR method are smaller than the symbols. 
Because the zero of free energy is arbitrary at each pressure, curves are vertically offset for clarity. Dotted lines are guides to 
the eye. c: Snapshots taken from umbrella sampling simulations biased to the stated packing fractions. Color scheme as in 
Fig. [3Ji. Scale bars = 50 nm. 



Remarkably, the resulting coexistence region encom- 
passes the grana packing fractions and compositions re- 
ported from many in vivo and in vitro experiments (e.g., 
filled and open squares in Fig. [5| . Other experiments re- 
port values of (p, <j>) within the model's fluid phases, but 
not far from the coexistence curve. Thus, our model of 
the LHCII-PSII protein system supports coexistence in 
many physiologically relevant conditions. 



III. DISCUSSION 

We have demonstrated that simulations of a sim- 
ple model of PSII and LHCII in stacked grana mem- 
branes, when configured to represent a generic grana 
membrane, can recapitulate many disparate and nontriv- 
ial experimental observations. These behaviors emerge 
spontaneously from the model's short-ranged interac- 
tions without the need for us to presuppose any par- 
ticular target assembled structures. Specifically, we ob- 
serve a distribution of C2S2, C2S2M, and C2S2M2 com- 
plexes (Fig. [2^,) ; LHCII-mediated intermembrane associ- 
ations between PSII supercomplexes (Fig. [3k— c) ; and co- 



occurence of disordered and crystalline-ordered regions in 
PSII- and LHCII-rich membranes (Fig. [3ji) . Importantly, 
all of these qualitative features appear to be common to 
grana membranes from many photosynthetic organisms 
and many growth conditions, although the quantitative 
details of course vary. To the best of our knowledge, this 
work marks the first reported computational investiga- 
tion of grana membranes to share these commonalities 
with experiment. 

Our principal prediction from numerical studies of this 
model is that the appearance of finite arrays of PSII 
and LHCII signals thermodynamic coexistence of disor- 
dered (fluid) and ordered (crystalline) phases. The phase 
boundaries we have computed further suggest that many 
physiological conditions lie at or near such coexistence. 
The experimental data reported in Ref. [28 , correlat- 
ing degree of array formation with the packing fraction 
and composition of grana membranes, offer the most con- 
crete opportunity for comparison. Membranes at con- 
ditions just inside our coexistence region, from plants 
grown in "ordinary" light, were found to exhibit a low de- 
gree (< 2%) of crystallinity. Membranes corresponding 
to conditions deep within our coexistence region, grown 
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FIG. 5: Phase diagram of crystal— fluid coexistence in 
the (<fi, p) plane. Each thermodynamic state {fit, p*) of 
phase equilibrium maps in this plane to two points and a cor- 
responding tie line. One point lies at the boundary of the ho- 
mogeneous fluid phase (low packing fraction and high LHCII 
content, shaded blue), the other at the boundary of the fully 
crystalline state (high packing fraction and low LHCII con- 
tent, shaded red). Example tie lines are drawn as dashed lines 
connecting these boundaries for selected values of p* . Coex- 
istence (shaded green) and pure-phase regions extend beyond 
the shaded areas determined by our limited data. In par- 
ticular, the coexistence region is expected to span the entire 
upper-right-hand quadrant. Filled and open symbols are low- 
light and ordinary-light phase points, respectively, calculated 
from data in Ref. [2S] (see SI Methods). 



in low light, showed significantly enhanced (22%) crys- 
tallinity. These consistent observations are indicated by 
the filled and open symbols in Fig. [5j The tie line we 
have computed suggests that a fully equilibrated system 
in ordinary light would possess greater crystallinity than 
observed in experiment, possibly reflecting high nucle- 
ation barriers and/or slow growth characteristic of dy- 
namics near coexistence (see Fig. S6). A recent study 
of low-light-acclimated membranes [29j reported arrays 
at protein packing fractions (p < 0.6) much lower than 
previously measured. Our calculations do not provide a 
direct way to resolve this apparent discrepancy, unless 
the different organisms feature substantially different en- 
ergy scales and/or modes of protein association. 

Experiments that systematically quantify dependence 
of crystallinity upon packing fraction at various protein 
compositions could further test or exploit our predic- 
tions. For example, grana membranes could be isolated 
from plants grown at different light intensities, generat- 
ing samples over a range of values of (j). Diluting these 
membranes with additional lipid [54] and measuring crys- 
tallinity via AFM or EM would allow construction of ex- 
perimental phase diagrams analogous to those we have 
computed. Matching experimental and theoretical re- 
sults in detail could determine appropriate values of the 



interaction parameters in our model. 

The well-understood phenomenology of phase transi- 
tions helps explain why the extent of PSII array forma- 
tion varies so dramatically between similar samples in 
experiment. The presence, extent, and number of arrays 
in the laboratory may be influenced by the characteris- 
tically slow dynamics associated with phase transitions. 
Nucleation of one phase from the other will be governed 
by rare structural fluctuations; crowding within grana 
membranes |34[I63| will likely make the subsequent repar- 
titioning of material between phases slow as well. More- 
over, grana isolation protocols that increase the protein 
packing fraction beyond the freezing transition densities 
for 2D hard discs (p « 0.7 [64 ) or hard rods {p 0.8 
[55]), such as BBY [53], could trap the isolated grana 
in jammed nonequilibrium configurations, further com- 
plicating experimental determinations of the equilibrium 
distribution of PSII arrays. Further computational work 
will be required to elucidate these dynamic effects. 

Proximity to phase coexistence could also contribute 
to substantial changes in thylakoid function observed to 
accompany modest changes in protein content and in- 
teractions. In addition to the low-light acclimation sce- 
nario discussed above, many regulatory processes asso- 
ciated with non-photochemical quenching, photoprotec- 
tion, and repair shift the system's position relative to 
phase boundaries in the p,(/>-plane: (1) State transi- 
tions involve transport of LHCII from PSII-rich to PSI- 
rich membranes [66) . reducing the local protein density 
and LHCILPSII ratio. (2) Photoinhibition-induced phos- 
phorylation decreases the diameter of grana stacks and 
breaks up PSII supercomplexes 67, 68 , changing the sys- 
tem size and composition. (3) The qE component of non- 
photochemical quenching may introduce an additional in- 
tramembrane attraction among LHCII [55] , affecting the 
relative stability of the fluid and crystal phases. These 
spatioregulatory processes are often interpreted as acting 
primarily over short length scales, tuning exciton fate by 
changing the relative distances between light harvesting 
sites, quenchers, and reaction centers. We suggest that, 
in addition, these processes may directly regulate global 
protein organization within the thylakoid membrane by 
inducing cooperative structural transitions. 

Suitably adapted, the model and computational frame- 
work developed in this study may help to clarify the 
mechanisms of such spatioregulatory changes. Direct 
analyses of grana-scale LHCII organization are difficult 
and rare, with the notable recent exception of Ref. [69] . 
With a few experimental parameters as input (namely 
<p and p or their equivalents), our methods can generate 
physiologically reasonable LHCII configurations in the 
presence of PSII in a physically grounded and unbiased 
fashion, and in sufficient quantity to enable statistical 
comparisons (e.g., Figs. S2 and S3). These ensembles of 
configurations could serve as the foundation for future 
studies of the many phenomena in which LHCII plays a 
key role. 



IV. MODEL AND METHODS 
A. Model geometry and energetics 

We represent protein complexes within membrane lay- 
ers of a grana stack as greatly simplified particles that 
can move in only two directions (x and y) and rotate 
only about the axis z perpendicular to these layers. Our 
simulations sample configurations of a pair of these two- 
dimensional systems, periodically replicated in x and y, 
coupled by stacking interactions between particles in dif- 
ferent layers. The model includes two particle species: 
isotropic disc-shaped particles represent trimeric LHC11 
complexes, and discorectangle-shaped particles represent 
PSII C2S2 supercomplexes. Particle shapes are simple 
approximations of the solved protein complex structures 
[3 [54], and particle sizes are assigned to be consistent 
with the protein structures and with previous coarse- 
grained models [33]. Specifically, each LHCII particle 
has a diameter ox = 6.5 nm, and each PSII particle has 
a rectangle width and cap diameter Dp = 12.0 nm and 
rectangle length Lp = 14.5 nm. The locations of the 
LHCIIs embedded in the PSII particles are modeled on 
Ref. [S3]. Full geometry details are given in Fig. SI. 

In addition to steric repulsions, our model includes 
three types of attractive interactions. The first two at- 
tractions (with energy scales esL-p and cml-p) act be- 
tween an interaction site on PSII and an LHCII particle. 
The interaction is square-well, with a distance cutoff of 
0.66ctl between the LHCII center and the binding site, 
or equivalently a distance cutoff between the LHCII edge 
and the binding site of 0.16ctl = 1-04 nm; this short- 
range distance cutoff ensures that each binding site can 
bind at most one LHCII at a time. This square well com- 
pletely describes the cml-p interaction between PSIIs and 
nonembedded LHCII particles. The attraction between 
PSIIs and embedded LHCIIs, esL-p, has an additional 
constraint: the angle difference between the long axes 
of the two PSII particles hosting the interaction sites 
must be < 30°. The depth of the square wells are set 
at esL-p = eml-p = 2 k^,T, except where noted. Tem- 
perature T fixed at 1 throughout, such that k^T is our 
reduced unit of energy. 

The other attraction involves two LHCII particles (free 
or embedded) in different layers. This interaction has 
an energetic minimum, with a value of cl-Lj when one 
disc completely eclipses another as viewed from above the 
layers. It vanishes continuously at r — ol, where r is the 
lateral distance between the centers of two LHCII discs. 
This dependence is plotted in Fig. [l]i and its functional 
form is given in the SI Methods. The energy scale is set 
to =4 fcp,T, except where noted. 

B. Monte Carlo simulations 

Using standard Metropolis Monte Carlo methods, we 
sampled equilibrium configurational distributions of sys- 



tems comprising Np = 2 x 128 PSII particles and a vari- 
able number Ax of free LHCII particles. Trial moves in- 
cluded translational displacements of a randomly chosen 
LHCII or PSII particle, and rotational displacements of 
a randomly chosen PSII particle. Translational moves 
were restricted to periodically replicated 2D surfaces, 
each with area A per layer, x- and y-components of 
translational moves were chosen from Gaussian distri- 
butions with zero mean and variances of 0.6 nm 2 for 
free LHCII displacements or 0.2 nm 2 for PSII displace- 
ments. PSII rotations were chosen uniformly on the in- 
terval (—10°, 10°). An "MC sweep" consisted on average 
of Ax LHCII displacements, ATp PSII displacements, and 
ATp PSII rotations. We performed simulations of closed 
systems (i.e., sampling the canonical ensemble) with sev- 
eral values of Ax and A, as given in Table SI. For each 
set of structural parameters, at least 7 simulations were 
initialized from independent fluid or partially crystalline 
configurations. 

To scrutinize phase behavior, we performed simula- 
tions in which A and Ax were allowed to fluctuate (i.e., 
sampling an osmotic ensemble |61j). These fluctuations 
were regulated by specified values of pressure p and chem- 
ical potential p,^, respectively. In addition to particle 
displacement trial moves, osmotic ensemble simulations 
included trial moves that change A (box area moves) 
and trial moves that change Ax (LHCII insertion and 
deletion moves). Acceptance probabilities for the latter 
moves were constructed to satisfy detailed balance, with 
the area per layer A scaled by the number of layers n; 
where necessary. 



C. Free energy calculations 

We determined free energy profiles F(p) by exploiting 
the fundamental relationship F(p) = — k#T In P(p) to 
the probability distribution P(p) for spontaneous pack- 
ing fraction fluctuations. Statistics of p were in turn 
determined by maximum likelihood estimation from a 
set of 139 osmotic ensemble simulations (£l cx .l = 0.1 
k B T relative to a standard state; p = 0.069, 0.070, 0.072 
k^T /nm?) with imposed harmonic bias potentials of the 
form Ubias(p) = k(p — po) 2 . Values of po and k were 
chosen to focus sampling on a series of small intervals 
spanning the range p = 0.625 to p = 0.852 (Fig. S5). Us- 
ing the MBAR method as implemented in the PyMBAR 
package [62] , these umbrella sampling results were pooled 
and reweighted to estimate F(p) at many pressures and 
chemical potentials. See SI Methods for details. 



D. Clustering algorithm 

PSII arrays were identified using a recursive clustering 
algorithm with three neighbor criteria, all of which must 
be satisfied for a PSII particle to be added to a growing 
array: (i) the angle difference between particle axes must 
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be < 15°; (ii) the component of the center-to-center sep- 
aration vector parallel to a particle axis must be < 26.5 
nm; (iii) the component of the center-to-center separa- 
tion vector perpendicular to a particle axis must be < 14 
nm. 
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